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^ ' Abstract 

The Local Fourier analysis (LEA) is a classic tool to prove convergence theorems 
for multigrid methods (MGMs). In particular, we are interested in optimality that 
is a convergence speed independent of the size of the involved matrices. For elliptic 
partial differential equations (PDEs), a well known optimality result requires that 
the sum of the orders of the grid transfer operators is not lower than the order 
^ ■ of the PDE to solve. Analogously, when dealing with MGMs for Toeplitz matrices 

in the literature an optimality condition on the position and on the order of the 

/ir zeros of the symbols of the grid transfer operators has been found. In this work we 

show that in the case of elliptic PDEs with constant coefficients, the two different 

[- — . approaches lead to an equivalent condition. We argue that the analysis for Toeplitz 

matrices is an algebraic generalization of the LFA, which allows to deal not only 
with differential problems but also for instance with integral problems. The equiv- 
alence of the two approaches gives the possibility of using grid transfer operators 
with different orders also for MGMs for Toeplitz matrices. We give also a class of 
grid transfer operators related to the B-spline's refinement equation and we study 
their geometric properties. This analysis suggests further links between wavelets 
and multigrid methods. A numerical experimentation confirms the correctness of 
the proposed analysis. 
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1 Introduction 



Multigrid methods (MGMs) are widely used for solving elliptic PDEs. The con- 
vergence analysis is usually done in the case of constant coefficients. Let us con- 
sider standard finite differences discretization for the following rf-dimensional 
problem 

d j2 9 

(- 1 ) 9 E -n^) = 9(x), xen = (o, if, q >i, 

i=i dV (1) 

boundary conditions on dQ, 

where x = (x\, . . . , xa). For simplicity of the presentation we assume the same 
meshsize h for each dimension, but generalization to h = (hi, . . . , ha) is possi- 
ble. Hence, approximating (1) by centered finite differences on a uniform grid 
of n subintervals of size h in each dimension, we obtain the linear system 

A n y = b (2) 

of order n d x n d . Neglecting the boundary conditions, the matrix A n is a rf-level 
Toeplitz matrix and it is banded at each level. 

Let Lh be a discretization of the differential operator in (1), then its Fourier 
transform is 

L{ U )=Y t li^ mw \ (3) 

for uj E \—-KJh,T\jh\ d , i 2 = —1, and 

L = -^— 7 / L(u) e-WM duj, (4) 

1 (27T) d J[-n/h,n/h]* K ' K ' 

where the operations between multi-indices are intended component wise and 
( • | • ) denotes the usual scalar product between vectors. Using centered finite 
differences of precision 2 and minimal bandwidth, the polynomial L(u) has 
degree 2q, that is the order of the PDE (1), and Lh is completely defined 
from its ci-dimensional stencil formed by the coefficients lj, with j t = —2q + 
1, . . . , 2q — 1 for i — 1, . . . , d. 

It is well-known that MGMs are optimal solvers for PDEs of the form (1), 
i.e., they require about a constant number of iterations varying h and each 
iteration has an arithmetic cost proportional to the matrix-vector product 
[20] . Such property is obtained imposing a well known condition on the order 
of the grid transfer operators: 

m r + m p > m, (5) 



where m = 2g is the order of the PDE, m r is the order of the restriction 
and m p is the order of the prolongation [3] , denoted as high frequencies order 
in [13,22]. The condition (5) follows from the LFA for the two grid method 
(TGM). In order to obtain more powerful grid transfer operators, that is in 
order to devise an optimal MGM, inequality (5) should to be satisfied strictly 
[20]. We note that the LFA does not consider the border effects, i.e., it assumes 
periodic boundary conditions or an infinite domain [3]. 

MGMs for multilevel positive definite Toeplitz matrices have been developed 
in the years looking only to the linear system (2), independently of the continu- 
ous problem [11,12,17,2,6,14]. A MGM for Toeplitz matrices was early defined 
in [11] using a powerful eigenvalue interlacing property with the matrices in 
the r algebra (the class of matrices diagonalized by discrete sine transforms 
of type I). This first proposal was extended to the multilevel case in [12,6,17]. 
Since Toeplitz matrices do not define an algebra and hence are difficult to 
manipulate, convergence results are provided using matrix algebra approxi- 
mations like r or circulant matrices having the same spectral distribution of 
the Toeplitz matrices. In other words, we require that the circulant or the r 
approximates share the same symbol of the original Toeplitz matrix [10,21]. 
In this paper we consider the circulant case. A MGM for circulant matrices 
was introduced in [19]. Furthermore in [2,1], generalizing the techniques used 
in [17], and using the Ruge and Stiiben theory [16] and the Perron- Frobenius 
theorem, a complete proof of the optimality of the V-cycle for multilevel circu- 
lant and r matrices was proposed. This analysis leads to a stronger condition 
with respect to the previous two grid analysis. 

In this paper we show that the techniques used in [11,12,17,19] represent a 
linear algebra generalization of the LFA [3,13], in the case of the Galerkin 
approach. Indeed, they lead to a condition analogous to (5), but on the or- 
der of the zeros of the generating functions of the grid transfer operators. 
The letter represents a wide generalization since the case of discretization 
of elliptic PDEs of order 1q corresponds to the case of generating function 
which are nonnegative (ellipticity), with a unique zero at zero of order 1q 
(consistency condition). In other words, the case of discretization of elliptic 
PDEs is a subcase of nonnegative symbols with unique zero at zero, which 
in turn represents the case when the algebraic problem is ill-conditioned in a 
subspace of low frequencies. Therefore by using the Toeplitz approach other 
cases can be considered including the case when the ill-conditioning arises in 
high-frequencies: we recall that the latter characterizes some integral problem 
related to signal/images restoration. We will show that considering the prob- 
lem (1), the LFA done in [3,13,22] and the Toeplitz approach [11,17] (which 
was introduced independently) are essentially equivalent. As already stressed, 
second approach is more general since it can be applied also when the zero of 
the symbol is not at the origin, or there exist several zeros (multiple sources of 
ill-conditioning). By the Galerkin approach, we have the only limitation that 



the restriction must be proportional to the transpose of the prolongation, but 
in this paper we will show that in practical implementations this condition is 
not necessary. More precisely, we will generalize the MGM for Toeplitz matri- 
ces to the case of a restriction different to the transpose of the prolongation. 
A first suggestion to consider the linear algebra tools for Toeplitz matrices as 
a generalization of the LFA for multigrid methods was given in [15,18]. 

In this paper we also define a class of grid transfer operators that satisfy the 
conditions in [13] but that are not interpolating. More in detail the considered 
operators are defined looking for the smallest support of the symbol for a fixed 
order and are related to the refinement equation of B-spline. We will give a 
geometrical interpretation of the operator of order 4 and it will be compared 
with the cubic interpolation. These B-spline grid transfer operators allow us to 
discuss some relations between wavelets and multigrid methods. Eventually, a 
numerical experimentation validates our proposals. 

The paper is organized as follows. Firstly in §2 we present the results in 
[13]. In §3 we describe the MGM defined in [11,12,17] for multilevel Toeplitz 
matrices using the zeros of the generating functions and we compare the two- 
grid analysis with the results in [13]. In §4 we generalize the TGM described 
in §3 to the case of restriction not necessary proportional to the transpose 
of the prolongation. In §5 we give a new class of grid transfer operators with 
minimum support for a fixed order and we show that such class is related to 
the B-spline. This allows, in §6, some observations about the relations between 
wavelet and multigrid methods. In §7 some numerical results validate the 
previous proposals both for Toeplitz non-differential problems and for PDEs 
with nonconstant coefficients. The final, §8 is devoted to some concluding 
remarks. 



2 The low and high frequencies order analysis 



We introduce a grid transfer operator that is not effective alone (it has order 
zero), but which is the basic tool for developing more powerful projectors. It is 
the classic down-sampling operator, called elementary restriction in [13] and 
cutting matrix in [11]. In the one-dimensional case, we set n = n^ > n^ > 
• • • > n® > 0, I e N, such that n (m) = (n {€) - (n {i) mod2))/2 and we define 
the down-sampling matrix K & G K" ' xn * as 



, lif 7=2Jfe-(nW + l)mod2, 

[*„(«>];,*=< k = l,...M l+1) - (6) 

otherwise, 



In the rf-dimensional case the down-sampling matrix is defined by tensor prod- 
uct as K ni i) = Km <g> K(i) ® • • • <g) Km. 



i 1 n 2 



Higher order grid transfer operators are defined by convolution with the down- 
sampling operator. The prolongation is P n (i){pi) = T n (i)(p i )K T (i) , while the 
restriction is R n (i)(ri) = K n ^)T n ^)(ri) that usually is the transpose, up to a 
constant factor, of the prolongation. The matrix T n (t) is the o?-level Toeplitz 
matrix generated by the function t. The Toeplitz matrices will be described 
in Section 3. The symbols r* and Pi should be trigonometric polynomials of 
low order to maintain the computational cost of the matrix vector product 
proportional to 0(N(n^')). 

Since more powerful grid transfer operators give a greater computational cost, 
it is important to find sufficient conditions such that we can decide for a fixed 
problem the cheapest grid transfer operators that allows to obtain an optimal 
MGM. This task can be done using the LFA [3,13] obtaining the condition 

(5). 

We define the set of all corners of x as Q(x) = {y\yj € {xj, n + Xj}, j = 
1, . . . , d}. With the change of variable x = uh, the set of all frequencies on 
the fine grid that correspond to the frequency u on the coarse grid is {z = 
y/h | y G Q(x)}. Moreover, according to the terminology in [11] we define the 
set of the "mirror" points of x as Ai(x) = Q(x) \ {x}. Since in the rest of the 
section we will consider only two grids, n will denote the fine grid. Moreover, 
for unifying the treatment, a generic grid transfer operator is denoted by 
B n (g) where g is multiplied by a factor 2 d when B n (g) is the prolongation, 
i.e., B n {g) = R n (g) or B n {g) = P n (2 d g). 

Definition 1 The Low Frequency order (LF) of a grid transfer operator 
B n {g) is the largest number s > for which 

g(x) = l + 0(\x\ s ), for \x\ — *• 0. 

Definition 2 The High Frequency order (HF) of a grid transfer operator 
B n (g) is the largest number s > for which 

g(y) = 0(\x\ s ), Vy <E M(x), for \x\ — > 0. 



For x = uh, \x\ — > means h — > since u is fixed. For h = (hi, . . . , hd) we 
can define |x| = max i= i v . .,d(|a;j|). 

For the grid transfer operators LF and HF are more general then classic in- 
terpolation order. 

Proposition 1 ([13])(%) If a restriction leaves all polynomials of degree s — 1 
invariant, then the LF of the operator is s. 



(ii) If a prolongation leaves all polynomial of degree s — 1 invariant, then both 
the LF and HF are at least s. 

For instance the linear interpolation has LF = HF = 2, while the cubic inter- 
polation has LF = HF = 4. 

Furthermore, we can derive the condition (5) from the following 

Proposition 2 ([13]) Given a constant- coefficient, linear differential opera- 
tor of order m, a necessary condition for non-increasing the high frequencies 
arising from a coarse grid correction with two grids it is 

lr + l P > m, (7) 

where 7 P and 7 r are the HF of the prolongation and of the restriction respec- 
tively. 

From Proposition 1 part (ii) the condition (7) is a generalization of the analo- 
gous condition on the interpolation order. The LF is important for the restric- 
tion thanks to Proposition 1 part (i), but it seems not necessary for the two 
grid analysis in Proposition 2. However, in [3] it is shown that for an efficient 
MGM a further condition is that both LF and HF are positive. This further 
request arises also from the Galerkin approach (see [22]) and is natural for 
obtaining an effective MGM. 

Eventually, we note that, since the grid transfer operation has to be computa- 
tionally cheap, the function g in Definitions 1 and 2 should be a trigonometric 
polynomial of low degree. Moreover, from Proposition 1 a good class of grid 
transfer operators should have at least LF> 0. Interpolating operators define 
a class with LF=HF. A further class of operators with a fixed HF and LF> 
will be described in Section 5. 



3 A MGM for Toeplitz matrices by generating functions 



In this section we briefly introduce the MGM defined and analyzed in [11,12,6,7,17,19,2,1] 
for the multidimensional r, circulant, Toeplitz and other matrix-algebras re- 
lated to trigonometric transforms. 

Toeplitz matrices arise from the discretization of convolution operators with 
a shift invariant kernel and hence not only from PDEs, but also from several 
other applications, e.g., image deblurring problems [2]. Toeplitz matrices are 
completely defined by the matrix size and the symbol also called generating 
function. Let / be a continuous function on M. d and having period 2-7T with 



respect to each variable, the Fourier coefficients of / are defined as 

a '-Wri„r f(x)e '' mdx - > e *- (8) 

Remark 1 With the change of variable x = uh, it holds a,j = lj and f(x) = 

From the coefficients {aj} one can build [21] the sequence {T n (f)} of multilevel 
Toeplitz matrices. Every matrix T n (f) is explicitly written as 

Uf)= E ••• E Hh,..,u)J^ ] ®---®J^ d] - 

|j'l|<n-l b'dl<"-l 

Here ® denotes the usual tensor product and Jjj^ G W nxn is the matrix whose 
entry (s, t) equals 1 if s — t — ji and is elsewhere, for i = 1, . . . , d. Many 
structural and spectral properties of T n (f) derive from its generating function 
/. Indeed, if / is real valued, then a_j = dj for every j and the matrices T n (f) 
are Hermitian for every n; if / is also non-negative but not identically zero 
then T n (f) is positive definite. 

Remark 2 The main difference between f(x) and L(u) is that uo denotes the 
frequency for the current discretization step h, information that seems to be 
lost in f , but that comes out from the matrix T n (f) regarding the current dis- 
cretization (h = l/(n + 1)). For instance, let Lh be the three point discretiza- 
tion of the Laplacian: then Iq = 2/h 2 and Z_i — l\ — —1/h 2 . On the other 
hand, in the algebraic approach for Toeplitz matrices the constant factor 1/h 2 
is moved to the right hand side (rhs) obtaining ao = 2 and o_i = ai = — 1. 
However the information of the order 2 of the Laplace operator is preserved 
since f(x) = 2 — 2cos(:r) vanishes at the origin with order 2. More in general, 
discretizing (1) with finite centered differences of minimal precision and mov- 
ing the coefficient l/h 2q to the rhs, by consistency, the symbol f(x) vanishes 
at the origin with order 2q. 

Convergence results for MGMs for PDEs are usually obtained neglecting the 
boundary conditions. In a similar way, MGMs for Toeplitz matrices are defined 
starting from matrix algebras like r or circulant. Imposing periodic boundary 
conditions in (1), the matrix A n in (2) is circulant. Circulant matrices are si- 

(n) 

multaneously diagonalized by the Fourier transform F n = -4= [e -y2/ > ]ij, where 
y(n) _ 2' K i/n, i — 0, . . . , n — 1. More precisely, the algebra of the circulant ma- 
trices can be formally defined as I A n \ A n = F n ■ diag(d) ■ F^ , deC n | where 

the vector d of the eigenvalues is equal to /(y^), y*-™*' = (Vo , • • • , 2/n-i)> an d 
a circulant matrix will be denoted by C n (f). In the <i-dimensional case the 
indices involved are multiindices, F n = F ni ® • • • <g) F nd has size N(n) and 
y*™) = y( rai ) x • • • x y( nd \ where x denotes the cartesian product. 



We do not consider boundary effects, thus we will discuss only the circulant 
case assuming periodic boundary conditions. In such case, in order to maintain 
the same circulant structure at each level, we have to start with n = n^ = 2 a , 
where a G N d . Moreover the grid transfer operators are defined as P n w (pi) = 
C n (i)ipi)K^ {i) and i? n (i)(n) = K„ w C nW (r;). In our case A n = C n (f) is singular 
since / vanishes at the origin which is a grid point. However, without losing 
generality, we assume A n nonsingular replacing / with its stabilized version 
that is by correcting A n by adding a special rank-one matrix. This correction 
is not consider here since it does not imply particular assumptions but it leads 
only to unnecessary complications in the notation [2]. 

Using the Galerkin approach, we must have R n (r) = P n {p) H , i.e., r = p, 
and A n /2 = P n {p) H A n P n {p) . Thanks to the structure of P n {p) we obtain that 
A n / 2 belongs again to the circulant algebra [19]. Thanks to the Ruge-Stuben 
theory [16], the TGM and the V-cycle convergence analysis can be split in 
two independent conditions, one on the smoother and the other on pi, for 
% = 0, ...,/ — 1, i.e., on the grid transfer operators. 

Remark 3 Several simple iterative methods, like relaxed Jacobi, satisfy the 
smoothing condition, therefore the main task is the study of the approximation 
condition for the grid transfer operators. 

In [19] the optimality of the TGM was proved for circulant matrices, under 
the following conditions on the grid transfer operators. 

Proposition 3 ([19]) Let the coefficient matrix be A n = C n (f) with f having 
a unique zero at x°. Defining P n {p) = C n (p)K^ and R n (r) = aP n (p) H , i.e. 
r = ap, a G M\{0} ; where p is a trigonometric polynomial non identically 
zero and such that for each x G [— n, 7i) d 



lim sup 

where 



p(yf 



fix) 



c < +oo, VyeM(x), (9a) 



E p(yf>^ (9b) 

y en(x) 
then defining A n /2 = aP n (p) H A n P n (p) the TGM is optimal. 

Proof. For a — 1 see [19]. For a ^ 1 it is enough to observe that the coarse 
grid correction CGC = I — P n (i? n A n P n ) _1 i? n A n is independent of a. □ 

In order to compare this result with the Proposition 2 we have to require 
p = 2 d r, thus the (7) becomes 2 / -/ r > m. We show the equivalence between two 
different convergence analysis for elliptic PDEs with constant coefficients: the 
LFA described in Section 2 and the analysis for Toeplitz matrices based on 
the zeros of the generating functions described here. 



Proposition 4 Let P n (p) = -Rn(2 d r) H 7 discretizing (1) by finite centered dif- 
ferences of order 2 and minimal bandwidth, the conditions (7) and (9a) are 
equivalent. 



Proof. By Definition 2 B n (g) has HF= s if and only if g(y) = with order s 
for all y G Ai(x). The discretization of an elliptic constant coefficient PDE of 
order m by finite centered differences of precision 2 and minimal bandwidth 
leads to A n = C n (f) (in the case of periodic BCs) with / vanishing at the 
origin with order m (see Remark 1). From condition (9a) p (or equivalently 
r) must be chosen such that p(y) = 0, for all y G ,M(0) with order 2 / -/ p > m. 
This is exactly the inequality in (7). □ 



The previous proposition shows that in the case of p = 1 d r and using the 
Galerkin approach, condition (9a) is a generalization of condition (7) to general 
problems not necessarily of differential type. The main difference between the 
two approaches relies in the coarse strategy. The results in [13] and summarized 
in Section 2 assume a discretization of the same PDE with the same formula 
at each grid. This imposes a right scaling of the grid transfer operators (i.e., 
B n (g) has LF > iff g(0) = 1). The latter is not necessary in the Galerkin 
approach adopted by the Toeplitz analysis, since the coarse matrix is defined 
as A n / 2 = aP n (p) H A n P n {p). Indeed the condition (9b) requires only p(0) ^ 0. 
More specifically, p can be defined up to a scaling factor since this gives only a 
different scaling of A n /2. However, the two approaches are comparable because 
from item 2 in Proposition 6, the coarse problem vanishes again at the origin 
and with the same order of the finer problem [17]. Using the PDE language, 
this means that for the Galerkin approach the linear system at the coarse grid 
is essentially (neglecting boundary conditions) the discretization of the same 
PDE with a formula of the same order. 



At the end of Section 2 we noted that Proposition 2 does not requires any con- 
dition on the LF of the grid transfer operators. The only interest on the LF 
could be deduced from Proposition 1, and mainly for the restriction. On the 
other hand, the TGM condition (9b) requires that the grid transfer operators 
have a positive LF (up to a scaling factor). This is exactly the same require- 
ment obtained in [22] for the Galerkin strategy and in [3] for an efficient MGM. 
In fact a condition LF= is equivalent to violate (9b) which implies p(x°) = 0. 
As a consequence, the associated grid transfer operators could fail to be full 
rank. The latter produces an increase of the ill-conditioning and could lead 
to singularity at the lower levels with a potential substantial change in the 
subspace related to small eigenvalues. 



MGM for Toeplitz matrices with a prolongation different from 
the transpose of the restriction 



In practical implementation the condition 

R n (r) = aP n (p) H (10) 

seems to be not necessary. The only request is that A n / 2 = R n (r)A n P n (p) 
is again positive definite for a recursive application of the algorithm. On the 
other hand, the condition (10) is very useful for a theoretical analysis, because 
if r 7^ p the coarse grid correction CGC — I n — P n (p)A~j 2 R n (r)A n is again 
a projector, but it is not longer unitary with respect to the scalar product 
< y, x>^ n = y H A n x, A n Hermitian positive definite, for all y, x e C n . For 
the well definiteness of a MGM, mainly to ensure that the same smoother is 
convergent also to the coarse levels, A n / 2 should be positive definite to apply 
recursively the algorithm. This condition is easy satisfied for p t ^ 0, r* ^ 
(not identically zero) and either both even or both odd, i — 0, ... ,1 — 1. More 
generally we could use p^i > with isolated zeros, i — 0, ... ,1 — 1. Therefore, 
the following generalization of Proposition 3 can be conjectured. 

TGM conditions. Let the coefficient matrix be A n = C n (f), with f having 
a unique zero at x° . Defining R n = K n C n (r) and P n = C n (p)K^ where p and 
r are trigonometric polynomials non identically zero and such that for each 

X E [— TV, it) 



\d 



lim sup 



r(y)p(y) 



/(*) 



c < +oo, Vy e M(x), (11a) 



where 



E r(y)p{y)^0. (lib) 

y&n{x) 

Then, defining A n / 2 = R n (r)A n P n (p), the TGM is optimal. 

These two conditions are motivated by the analysis in the previous section 
and by the following Proposition 5 that extends Proposition 4 to the case of 
r ^ ap. Moreover, the numerical experiments in Section 7 will validate these 
conditions. 

Proposition 5 Discretizing (1) by finite centered differences of order 2 and 
minimal bandwidth, the conditions (7) and (11a) are equivalent. 

Proof. The proof is analogous to that of Proposition 4. It is enough to observe 
that if r and p vanish at y with order 7 r and 7 P respectively, then rp vanishes 
at y with order 7 r + 7 P . □ 

We provide a further result useful to implement the corresponding MGM. 
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Proposition 6 Let A n = C n (f), P n (p) = C n (p)K^, R n (r) = K n C n (r), with 
f,p,r trigonometric polynomials andp,r satisfying conditions (11). Then 

(1) A n / 2 = R n {j-)A n P n {p) coincides with C„/ 2 (/) where 



f(x) = 7^ E r (y)f(y)p(y)> xe[-n,ir) d . 



2 d 



:i2) 



y€Cl(x/2) 



(2) If x° G [—7i,ii) d is a zero of f , then y° = 2a;°mod27r is a zero of f. 
Moreover the order of the zero y° of f is exactly the same as the one of 
the zero x° of f. 

Proof. The essentials of the proof in the case of r = p can be found in [19]. For 
r^pwe can proceed similarly. We sketch the main steps for the one dimen- 
sional case and, at the end, we extend it to the multidimensional case, mainly 
for emphasizing the algebraic interpretation of the frequencies packaging used 
in the LFA. 



The main relationship is 



K F 



1 

V2 



H n/2 | f n j2 



that implies 



A n/2 = K n C n (r)C n (f)C n {p)K : 



1 r 

2 



F n /2 I F n / 2 



(li %=0,...,n-l(r/p(-^) 



n 



b n/2 
b nl2 



*-*• '-**. 



rfp 



n/2 



TV 



n/2- 



that is the (12) for d = 1. 



In the multidimensional case K n F n = 2 d//2 -F„/ 2 G„, where F n = ®j = i F n ./ 2 
and G n = ® d j=l ([1 1] ® 7 n . /2 ). Therefore G n D n (rfp)G T n = D n/2 (2 d f), where 
D k (h) = diag ^ j ^ k _ e (h(2nj/n)) and / is defined by (12). 



The claim in item 2 is a consequence of item 1 and of relations (11). 



□ 



Thanks to Proposition 5, the TGM conditions in (11) give a complete gen- 
eralization of Proposition 2, also for non-differential problems since the case 
of generating functions vanishing at points different from the origin is also 
included. Therefore the analysis based on the zeros of the generating func- 
tion can be considered an algebraic generalization of the LFA also to non- 



11 



differential problems. For instance, some discretized integral problems have a 
generating function vanishing at ne, e = (1, . . . , 1) T G M. d , or more generally 
at some nx, x t G {0, 1}, i — 1, . . . , d, ||a;||oo = 1, with order 2q. In this case 
Hg(x) = 2~ dq n1j =1 (l — e~ 1Xj )\-2i e^'h-l satisfies the conditions (11) and there- 
fore it defines an optimal TGM. We note that R n = K n T n (ii q ) is an high-pass 
filter and then it projects into the high frequencies. However, thanks to Propo- 
sition 6, the zero at the next level moves to the origin and at the coarser grids 
the problem becomes spectrally equivalent to the discretization of a constant 
coefficient elliptic PDE. 

Finally, we recall the V-cycle optimality conditions for Toeplitz matrices given 
in [2] for d = 1 and in [1] for d > 1. 

Proposition 7 ([2,1]) Let A n (i) = C n (i)(fi) be the coefficient matrix at the 
level i, for i = 0, . . . , /, with fi having a unique zero at x®. Defining A n (.+i) = 
R n (i)(ri)A n (i)P n (i)(pi), P n (i)(pi) = C n (i){pi)K^ and R n (i)(ri) = aP n ( l) (p i ) H , i.e. 
Ti = api, a G M\{0} 7 where Pi is a trigonometric polynomial non identically 
zero and such that for each x G [—it, n) d 



lim sup 



Pi(y) 



fi(x) 



c < +oo, \/yeM(x), (13a) 



where 

E Mv) 2 >Q- (isb) 

yan(x) 
Then the V-cycle is optimal. 

We observe that (13a) defines a stronger condition on the order of the grid 
transfer operators with respect to the condition (9a). On the other hand, 
choosing c = in (9a) or in (11a), which is equivalent to require that the (7) is 
satisfied strictly, is usually enough to obtain an optimal K-cycle as numerically 
shown in [19,20,2]. Following the same analysis done for Proposition (3), the 
Proposition (7) could be generalized to the case of r ^ p and also applied 
to non-constant coefficients PDEs. In this last case, condition (13a) could be 
rewrite as 7 r + 7 P > 2m. 



5 B-spline grid transfer operators 



From the discussion at the end of Section 2, the HF is more important than 
LF: in fact, for the latter it is enough require LF> 0. A class of grid transfer 
operators having HF= m and LF> can be defined by 

^(x) = U\^ L 2—) ■ ( 14 ) 
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Therefore every grid transfer operator with HF= m has a generating function 
of the form (p m (x)v m (x) such that v m (x) ^ for all j6M(0) and v m (0) = 1. 
With v m = 1 we obtain a class of projectors with minimal support for a fixed 
order m. We note that v2y? m is the symbol of the B-spline of order m in the 
multiresolution analysis (MRA) [8]. This technique is based on hierarchies of 
nested spaces Vj C Vj +i , j G Z, defined through a basis generated (for d — 1) 
by translations and dilations f3(2 : >x—k), k G Z, of a single scaling function /3. A 
scaling function satisfies an equation of the type /3(x) = V2J2kez hkP(2x — k), 
which expresses the nestedness of the spaces Vj = span{/5(2- J x — k),k G Z}. 
Let /3(v) be the Fourier transform of f3(x). Then 

(3(v) = H(v/2)P(v/2), (15) 

with H(v) = l/V2Y.kezhk e ^ lkv - The function H{v) is called the symbol of (3. 
A common application of the MRA is to approximate a high resolution / G Vj 
by a coarser function / G 14 with k < j without losing a lot of information. 

In the one dimensional case a simple scaling function is the Haar-function 
/3(x) = 1 for x G [0, 1) and zero otherwise, which satisfies the refinement 
equation Si(x) = S\(2x) + S\(2x — 1). The Haar-function Si is the simplest 
B-spline of order m — 1, and {Si(x — k) : k G Z} is an orthonormal basis of 
L 2 (1R). The B-spline of order m can be defined by S m — Si * S m -i, where * 
is the convolution operator. For instance S2 = S± * Si is such that ^(x) = 
l/2(S 2 (2x) + 2S 2 (2x - 1) + S 2 {2x - 2)) and its translated S 2 (x) := S 2 {x + 1), 
S 2 (x) = 1 — |x| for x G [—1, 1] and zero otherwise, is known as the hat function. 
We remark that S 2 = 2ip 2 and more generally S m = 2ip m . 

The functions S m are not centered, but they can be easy centered as previously 
done for S 2 : instead of v m = 1 we take the shift v m = e 1:r l-TJ. In this way we 
define a class of centered projectors such that the symbol of order m is 



s) = n : eixlfl - ( l6 ) 




The (j) m have HF= m and LF= 2. As previously observed, we note that a 
good class of grid transfer operators should have a high HF, while for the LF 
the only request is LF > 0. For d — 1, the m , can be obtained using the 
Tartaglia's triangle as in Table 1. For d = 2 we take the tensor product of the 
one dimensional stencil and so on for d > 2. 

The class of projectors defined in (16) is a scaled generalization of p 2 k{x) = 
V / 2(l + cos(x)) fc proposed in [2] for even functions (zeros of order m = 2k) 
in the one dimensional case, indeed p 2 k = 2 1 / 2 ~ fe 2 fc- However, scaling factors 
does not change the effectiveness of the projector for the Galerkin approach. 
We remark that if m is odd than the grid transfer operators related to <p m are 
not symmetric for vertex centered discretization, while they are symmetric for 
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Table 1 

The refinement coefficients hk 7^ 0, k £ Z for 2 r ' 



in the one dimensional case. 



m 


fc-2 


h-i 


/i 


hi 


h 2 


1 




1 


1 






2 




1 


2 


1 




3 


1 


3 


3 


1 




4 


1 


4 


6 


4 


1 



cell centered discretization [22]. For instance [13 3 1] is the linear interpolation 
for cell centered discretization. 

We consider for simplicity the one dimensional case, but the following obser- 
vations are true also for d > 1. For vertex centered discretizations, we are 
interested in tf2k- It is easy to prove that P n {4>2) is the linear interpolation. 
Which is the geometrical meaning of P n (04)? Which are the relations between 
Pn{4>4) and P n (g c )l Answers to these questions will be given in the next sub- 
section. 



5.1 A quadratic prolongation of order 4 



In this subsection we give a geometric interpretation of P n 
HF= 4 like the cubic interpolation. 

The simplest but useless prolongation is 



which has 



(A) 



p(x) 



1. 



Without losing in generality we consider n odd. For y G C n and x G C - ^ - , 
y = P n (l)x = K n x does not reconstruct constant functions not identically 
zero. Particularly, the choice (A) does not provide a good approximation for 
the odd components. Therefore, in the standard MGM, P n (p) is frequently 
chosen as the linear interpolation 



(B) 



p(x) 



1 + cos(x) 



Remark 4 The choice (B) compared to the choice (A) leaves unchanged the 
even components but reinforces the odd components, which are not well ap- 
proximated by the choice (A), with a linear interpolation. 

When the choice (B) is ineffective, it is usually replaced with the cubic inter- 
polation. An alternative is given by P n (<p4) which has HF= 4 like the cubic in- 
terpolation but a smaller support. This prolongation follows a strategy similar 
to that used for deriving choice (B) from (A): it leaves the linear interpolation 
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for the odd components and reinforces the even components. From (16) 

d 

4 (a;)=4- d n( 1 + cos (^)) 2 - 
i=i 

Therefore, in the one dimensional case, r(x) = (1 + cos(x)) 2 /4 and 

(C) p(x) = (1 + cos(x)) 2 /2. 

Remark 5 With respect to the choice (B), this choice leaves unchanged the 
odd components but reinforces the even ones with a quadratic approximation: 

f(x fc + x fc+ i)/2, J = 2k + 1, 

y j = I fc = l,...,ra, (17) 

[ (x fc _i + 6x fc + x fc+ i)/8, j = 2k, 

where we assume xo = x n+ i = 0. 

The approximation for the even components of y is obtained taking the mid- 
dle value of a quadratic rational Bezier curve defined from the three points 
{xfc_!,Xfc,Xfc + i}. The Bernstein polynomial of order n is defined as 

Bt\t) = 
A quadratic rational Bezier curve has the expression 




C(t) 



n*<*iB?\t) 



where bj are the control points and u>i are the associated weights for i — 0,1, 2. 
Let bj = x fe+ j_! for i — 0, 1,2, co 1 = 3/2 and u = u 2 = 1/2, then C(|) = 
(xfc_i + 6x fc + x fc+1 )/8 which is the same of (17). In Figure 1 the previous 
quadratic approximation is shown for the computation of y^ with j = Ik. 
Furthermore in Figure 2 we compare the values obtained in the finer grid 
with the choice (B) (linear interpolation) and with the choice (C) (quadratic 
approximation) . 

We note a different philosophy between the cubic interpolation and this quadratic 
approximation. The cubic interpolation reconstructs exactly the polynomial 
of degree at most three, while the quadratic approximation does not. However, 
for the even nodes the cubic interpolation preserves the exact value as in the 
coarse grid. This can be useful for the TGM or when the coarse solution is 
well approximated. But when this is not the case, like for the V^-cycle in some 
applications, it could be better to take an approximation that is an average 
at the neighboring nodes. The underling idea is that it is not useful to take 
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x* 




Fig. 1. Even components in the finer grid computed with the choice (C) 




Fig. 2. Computation of the points y« for i = 1, . . . , 5 in the finer grid using the linear 
interpolation (line) and the quadratic approximation (dotted). 

a powerful prolongation and a poor restriction, because in this case we inter- 
polate the solution of a coarse problem that does not represent well the finer 
problem. 



The generating function of the cubic interpolation is g c 



v(x) = 2 — cos(a;) 
the one of B„ (6a) 



[x] = (p4.(x)v(x), with 
9 0-1] and 



The stencil of -B n (g c ) is 32 [—1 9 16 
is 4r[l 4 6 4 1]. In addition to the different philosophy 
previously emphasized, we observe the following mainly differences. 

Remark 6 B n (g c ) has HF=LF=4, while B u (6a) has HF=4 and LF=2. 

Remark 7 From a computational point of view, the two stencils have the 
same number of nonzero elements. Hence they have the same computational 
cost for the projection of a vector between coarse to fine or fine to coarse grids. 
The main difference is that the stencil of the cubic interpolation has a lager 
support. This implies that, mainly for d > 1, B n (6^) defines a MGM that is 
computationally more efficient at the coarser grids. Indeed, using the Galerkin 
approach, the coarse matrix has a lower bandwidth. Moreover, _B„(0 4 ) requires 



16 



less boundary points, property that is very useful for a parallel implementation. 
In this way less communications are required among the nodes and we can 
employ a colored Gauss-Seidel with a smaller number of colors, increasing the 
parallelism degree. 

From the previous remarks, as it will be confirmed by the numerical experi- 
mentation in Section 7, the cubic interpolation usually converges within less 
iterations with respect to the choice (C) (see Remark 6), but both have the 
same asymptotic behavior since they have the same HF. Therefore, thanks to 
Remark 7, i? n (0 4 ) could be a good alternative to B n (g c ), mainly for parallel 
implementations or V-cycle MGMs. 



6 A comparison between MGMs and wavelets methods 



In this subsection we recall some approximation properties of the m defined 
in (16), or equivalently of the f m in (14). Moreover we discuss some relations 
between wavelets and multigrid methods. 

Unfortunately {S m (x — k) : k G Z, m > 1} is not an orthogonal system. 
However, the S m satisfy the quasi- interpolant Strang-Fix condition [8]: 

f^j S m (2kn) = 0, keZ\{0}, s = 0,...,m-l. (18) 

It follows that the polynomials of degree at most m — 1 are contained in the 
space Vq = span{S' m (a; — k) : k G Z}. The same property is also usually 
expressed in terms of vanishing moments. Let ifj m be the wavelet associated 
to S m . The first m — 1 moments of ip m vanish, i.e. 

x s ip m (x) = 0, s — l,...,m — l. (19) 

-oo 

Starting from the orthogonality condition 

\H(v)\ 2 + \H(v + n)\ 2 = l, (20) 

where H(v) is the symbol defined in (15), and imposing the vanishing of the 
moments, Daubechies defined orthogonal wavelets [9] . The moment condition 
(19) is equivalent to require H(x) = (p m (x)v(x), such that v(x) is a trigono- 
metric polynomial with u(0) = 1. This means that orthogonal wavelets are 
obtained imposing orthogonality to S m . For instance, for m = 2we obtain the 
Daubechies wavelet of order 2 with scaling coefficients h = -K= [1 + y/3, 3 + 

y/3, 3 — a/3, 1 - a/3] and symbol (p 2 (x)(l + \/3 + (1 - \/%)) e~ [x )/2. 
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Considering transfer grid operators for multigrid methods, a question now 
seems to be natural: "Is the orthogonality condition (20) necessary or the 
moment condition (19) is sufficient?" The answer is that the orthogonality 
condition (20) is not necessary. Some motivations were firstly given in [4]. 
Indeed, in the MRA we would have Vo = Vi© Wi, while for multigrid methods 
we would have Q = Range{P} © NullSpace{P} where Q h is the fine grid [5]. 
The error e in Q h can be decomposed as e = s + t, with s G RangejP} and t e 
NullSpace{P T } (we consider the Galerkin condition R = P T ). We note that 
CGCs = 0, while CGCt = t. Therefore, if after the pre-smoothing step t — 0, 
then the TGM converges in one iteration. Clearly this is a too strong condition. 
However Range{P} is spanned by smooth functions while NullSpace{P T } is 
spanned by oscillating functions, thus if we have a "good" smoother t ~ 0. 
Consequently, to define a good grid transfer operator we are interested only 
to the scaling function and we do not use the wavelet since the smoother has 
already reduced the error in the high frequencies. Eventually, the multigrid 
methods is an iterative method, hence it is not necessary to have convergence 
in one iteration but the solution can be substantially improved iterating. These 
considerations show that for multigrid methods the orthogonality condition 
(20) is not necessary. The only necessary condition concerns moments in (19), 
that is equivalent to the HF and to the factorization (p m v m . 

Concluding, for Toeplitz linear system relations between wavelets and multi- 
grid methods was already investigated in [7]. In such paper the authors proved 
the TGM optimality using Cohen, Daubechies and Feauveau (CDF) 9/7 biortho^ 
onal wavelets for generating functions having a zero of order 4 in the origin. 
Thanks to the previous comments this is expected since CDF 9/7 have both 
moments of order 4. However the interesting fact is that the proof in [7] is done 
directly in the Toeplitz class and not in the algebra case. Moreover, comput- 
ing a class of compactly supported biorthogonal wavelet systems GCDF with 
specified vanishing moments for scaling functions, the authors show numeri- 
cally that for moments of order k the TGM is optimal for generating functions 
having a zero of order at most Ik. This is exactly the condition (9a) since the 
symbols have the form H(x) = tp(x)v(x). 



7 Numerical Experiments 



For the sake of simplicity, we consider some tests in the ID case. The re- 
sults in the multidimensional case are very similar. The numerical experiments 
are done using Matlab 7.0. We fix the following parameters. The MGMs are 
stopped when the relative l-i norm of residual is lower than 10~ 9 . The coarse 
matrices are defined using the Galerkin approach. The coarser problem has di- 
mensions 7x7. These means that A n a) = T n w (/,) + low rank, for i — 1, . . . , /, 
where the f\ are defined according to (12). The right hand side is obtained 



from the exact solution Xj = j/n, j = 1, . . . ,n. The pre-smoother and the 
post-smoother are one step of relaxed Richardson with relaxation parameter 
equal to 1.5/||j^|| 00 and 1/H/iH^ according to [1]. 

In Section 4, Proposition 5 gives a validation of the TGM conditions (11) 
for PDEs. Here, in the first test problem, we give a numerical validation of 
the conditions (11) in the case of discretized integral problems. We show that 
sometimes the choice r ^ p can be very useful also for the MGM for Toeplitz 
matrices. In the second test problem, we show that the V-cycle optimality re- 
sult for Toeplitz matrices in Proposition 7 can be applied also to non-constant 
coefficients PDEs. 



7.1 An integral problem 



We consider the discretization by the rectangle quadrature formula of the 
following Fredholm operator of first kind: 



g(x) = [ t(x - 9) f(9) d6, x, 6 e 



*j 



where / is the input object, t is the integral kernel of the operator, also called 
point spread function (PSF) and g is the observed object. In the discrete case, 
when zero Dirichlet BCs are used and the PSF is shift invariant, the above 
approximation gives rise to the system A n f = g, where A n = T n (z) with 
z(tt) = and positive elsewhere. According to the analysis in Section 4, the 
generating function of a grid transfer operator of order s will be /i s at the finer 
level and (f) s at the lower level. We denote by 5 r and S p the order of the zero 
of the restriction and of the prolongation, respectively. 

First of all, we give a numerical evidence of the relevance of conditions (11). 
We consider z(x) = (2 + 2cos(x)) 3 , which has a zero in it of order 6. From 
Proposition 3 S r = S p = 2 is not enough for an optimal TGM, while it is 
necessary to set S r = S p = 4. However, if we allow S r ^ S p then, from conditions 
(11), the optimality of the TGM is guaranteed with 5 r = 2 and 5 P = 4. This 
is confirmed by the numerical results in Table 2. In such table we report the 
number of iterations required by the TGM for converging for different orders 
of the grid transfer operators, when increasing the problem size. 

In real applications usually more than two-grids are used. The TGM optimality 
conditions are not enough for the V-cycle, but they give good estimations for 
the H^-cycle [20]. Defining z(x) = (2 + 2 cos(a;)) 2 , which has a zero in 7r of order 
4, also the choice S r = 5 P = 2 gives an optimal TGM. However, in Table 3 we 
see that for the H^-cycle there is a large reduction in the iteration number for 
S r = 2 and S p = 4, with respect to 5 r = 5 P = 2. Furthermore using 5 r = S p = 4 
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Table 2 

TGM iteration numbers varying the order of the grid transfer operators and the 
problem size n for the integral problem z(x) = (2 + 2cos(x)) 3 . 



n 


S r = 2 


5 r = 2 


S r = 4 




5 P = 2 


<5 P = 4 


5 P = 4 


15 


219 


65 


51 


31 


607 


72 


52 


63 


1501 


76 


51 


127 


> 2000 


77 


50 


255 


> 2000 


78 


49 



implies a negligible reduction of the iteration number, with respect to the 
choice 5 r = 2 and S p = 4. 

From a computational point of view, we should investigate the structure of the 
coefficient matrices at each level for the previous choices of the grid transfer 
operators. Indeed, at each level the main computational cost is related to 
the matrix vector product with the matrices A n (i), for i = 0,...,/. In the 
last example with z(x) = (2 + 2cos(x)) 2 , we use the following grid transfer 
operators of order 2s: r = Po = (2 — 2cos(:r)) s and r i = , p i = (2 + 2cos(a;)) s for 
i — 1, ... ,1 — 1. For S r = S„ = 2 we have A„u) = Tj^(z), for i — 1, . . . , I, where 



z(x) 



(2 — 2cos(x)) . For 5 r = 2 and 5 P = 4 we have A r 



(i) 



2 i T (i) ( 



z) + 



Cit\e{ + Qe„e„, where e, is the jth vector of the canonical basis and q = 
[T^(ri)A n (»)T^(pj)] 2 ,2 — 2 z -6. Therefore, the matrix vector product has about 
the same computational cost for both choices. On the other hand, for S r = 
5 P = 4 the matrix A n w is Toeplitz plus a 4 rank correction and moreover the 
bandwidth of the Toeplitz part is not longer 5, but it becomes 7. Obviously, 
this fact increases the complexity and the computational cost of the matrix 
vector product. Moreover, for using the Gauss-Seidel smoother with a coloring 
strategy we need 4 colors instead of 3 colors as for T n (z), losing a degree of 
parallelism. The previous considerations are enhanced in the multidimensional 
case. Indeed in the two dimensional case the bandwidth of each block moves 
from 5 to 7, and also the block bandwidth moves from 5 to 7. 



For preserving the Toeplitz structure at each level, a different cutting matrix 
proposed in [2] can be used. The main idea is in the changing of the coarse 
problems size in order to neglect in some way the boundary effects that give 
the low rank corrections. However, in this way we lose some information and 
the iteration number slightly increases even if the general (optimal) behavior 
is preserved in the numerical experimentations. 
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Table 3 

VF-cycle iteration numbers varying the order of the grid transfer operators and the 
problem size n for the integral problem z(x) = (2 + 2cos(x)) 2 . 



n 


S r = 2 


S r = 2 


S r = 4 




5 P = 2 


5 P = 4 


5 P = 4 


31 


25 


23 


22 


63 


32 


23 


21 


127 


35 


23 


21 


255 


37 


23 


20 


511 


37 


23 


20 



7.2 A differential equation with nonconstant coefficients 



We consider the following equation 




g{x) 



xe(0,l), 



(21) 



with nonconstant a(x) and order m = 4. 

In this subsection we consider the V-cycle, that is cheaper than the W^-cycle 
in vector and parallel implementations. Therefore we need a more powerful 
smoother and we replace the weighted Richardson with Gauss-Seidel. The 
generating function of the grid transfer operator is the same for each coarse 
problem. We test several combinations of <pk, k = 2,4 and the cubic interpo- 
lation g c . 

We note that in practical implementations it is usually required that (7) is 
satisfied strictly. This condition in terms of generating functions is equivalent 
to require that c = in (11a): the latter is a generalization of a similar 
condition already numerically observed in [2] for r = p and when considering 
the l/-cycle. Indeed, from Table 4, we see that the choice (02, 02) for the 
couple (restriction, prolongation) is not effective. We recall that 02 is the 
linear interpolation. For obtaining an effective V-cycle it is enough to increase 
only the order of the prolongation (i.e., 7 r = 2 and 7 p = 4). However, more 
stable results can be obtained for 7 r = 4 and 7 P = 4 according to Proposition 
7. 

Eventually, we compare the two prolongations generated by 04 and g c . In terms 
of number of iterations the function g c has to be preferred with respect to 04 
according to Remark 6. On the other hand, 04 has a 5-point stencil while g c has 
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Table 4 

y-cycle iteration numbers varying problem size n and a(x) for the differential 
problem (21) (g c = cubic interpolation). 



restriction 


02 


02 


02 


04 


04 


02 


02 


02 


04 04 


prolongation 


02 


04 


9c 


04 


9c 


02 


04 


9c 


04 9c 


n 




a(x) = 


C x ' 




a 


>) = 


(x - 


0.5) 2 


15 


14 


9 


9 


7 


7 


15 


10 


10 


9 9 


31 


32 


11 


13 


10 


9 


33 


13 


17 


10 11 


63 


60 


17 


15 


14 


9 


61 


17 


24 


13 11 


127 


98 


27 


20 


18 


12 


101 


26 


27 


17 13 


255 


151 


38 


27 


22 


16 


155 


35 


29 


20 16 


511 


215 


48 


34 


26 


20 


221 


44 


36 


24 19 


1023 


276 


57 


44 


29 


22 


284 


53 


46 


27 22 



a 7-point stencil. This means that the choice (02, g c ) leads to coarse matrices 
having a bandwidth equal to 7 like the choice (04,04), while (04, g c ) leads to 
a bandwidth equal to 9. Therefore, the right comparison should be between 
(02,,9e) and (04,04). From Table 4, it is evident that the second one has to 
be preferred according to Remark 7. Moreover, considering the discretization 
near the boundary, 0& requires less boundary points than g c . 



8 Conclusions 



Considering elliptic PDEs with constant coefficients, we have shown the equiv- 
alence between the LFA and the analysis based on the zeros of the generating 
functions of Toeplitz matrices. This equivalence has two implications. The 
first one is that the techniques used for Toeplitz and Circulant matrices allow 
to extend the LFA also to non-differential problems (e.g., integral problems). 
The second one is that it suggests to choose the restriction different from the 
prolongation also for MGMs for Toeplitz linear systems. The generalization of 
the MGM for Toeplitz matrices proposed in Section 4 replaces the Galerkin 
conditions with the following: 

(1) A n/2 = R n {r)A n P n {p) 

(2) r not necessary equal to p, but such that p > 0, r > and both even or 
odd (such that A n / 2 is again positive definite). 

We have given a class of grid transfer operator with minimal support for a 
fixed HF and the geometrical interpretation of the operator with HF = 4. Such 
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class is related to the symbol of B-spline, giving the possibility to discuss some 
useful and suggestive relations between wavelets and multigrid methods. 
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